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ABSTRACT 


We  investigate  the  ISAR  imaging  model  using  two  different  approaches  the 
Fourier  and  the  wavelet  transform.  Starting  from  the  weak-scatter  far-held  model  we 
explain  why  this  approximation  it’s  a  Fourier  Transform.  We  investigate  the  scatter¬ 
ing  mechanisms  from  different  surfaces.  We  also  analyze  and  derived  the  expressions 
of  the  scatter  from  sphere.  We  briefly  discuss  the  radon  transform  and  why  it’s  not  a 
good  approach  to  our  problem.  We  investigate  the  results  using  this  time-frequency 
methods,  creating  images  using  synthetic  data. 
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I.  INTRODUCTION 

A.  SYNTHETIC  APERTURE  RADAR  IMAGING 

Radar  is  an  electromagnetic  instrument  used  for  detection  and  localization  of 
targets,  such  as  aircraft,  ships  and  ground  vehicles.  It  transmits  energy  and  receives 
the  backscattered  held  energy  from  the  target  [Ref.  4],  The  received  signal  is  then 
processed  to  extract  information  about  bearing  and  distance.  The  radar  as  a  ra¬ 
dio  frequency  sensor  can  perform  at  long  ranges  with  high  accuracy  independent  of 
weather  conditions,  in  day  or  night. 

With  the  advent  of  powerful  digital  signal  processing  algorithms,  multidimensional 
signal  analysis  and  inversion  in  imaging  systems  may  now  be  formulated  via  more 
concrete  theoretical  principals.  These  principals  will  be  the  central  topic  of  this  the¬ 
sis:  radar  technology  and  the  powerful  tools  of  digital  signal  processing  (DSP)  can 
extract  more  information  from  the  received  signal  than  just  range  and  bearing.  In 
our  case,  we  want  to  be  able  to  produce  images  from  moving  targets. 

The  original  radar  systems  measured  range  to  a  target  (radar  scatterer)  via  time  de¬ 
lay,  and  direction  of  a  target  via  antenna  directivity.  In  the  past,  Doppler  shifts  were 
used  to  measure  target  speed  only,  until  it  was  discovered  that  Doppler  shifts  could  be 
processed  to  obtain  fine  resolution  in  a  direction  perpendicular  to  the  range  or  beam 
direction.  Through  this  latter  concept,  often  credited  to  Carl  Wiley  of  Goodyear 
Aerospace  in  1951,  it  was  found  that  two-dimentional  images  could  be  made  of  the 
targets  and  of  the  earth’s  surface  using  radar.  The  method  was  termed  Synthetic 
Aperture  Radar  (SAR),  referring  to  the  concept  of  creating  the  effect  of  a  very  long 
antenna  by  signal  analysis  [Ref.  2].  The  ability  to  view  or  capture  a  scene  improves 
with  a  larger  aperture  (in  a  binocular  or  camera),  a  larger  radar  antenna  aperture, 
or  larger  acoustic  transducer  aperture;  the  key  to  better  vision  is  a  larger  aperture. 
Unfortunately,  it  is  extremely  difficult  to  either  build  or  maintain  a  physically  large 
aperture  radar  system. 
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The  principal  idea  behind  SAR  is  to  synthesize  the  effect  of  a  large  physical  aperture 
radar,  whose  construction  is  infeasible.  The  significance  of  a  larger  aperture  becomes 
evident  in  the  following  example: 

The  lateral  or  cross-range  resolution  of  a  D  =  1  meter  diameter  radar  antenna  with 
wavelength  A  =  1  meter  at  the  range  R  —  1000  meters  is 


RX 

Lateral-resolution  =  =  1000  meters, 


(1.1) 


which  is  very  poor.  Yet,  based  on  the  SAR  Theory  and  signal  processing,  if  we  move 
this  small  1  meter  aperture  radar  along  a  synthetic  aperture  with  length  De//=1000- 
meters  {Deff  stands  for  effective  radar  diameter  or  aperture),  then  the  lateral  reso¬ 
lution  with  wavelength  A  =  1  meter  at  the  range  R  =  1000  meters  becomes 


Lateral-resolution  = 


RX 


D, 


eff 


—  1  meter, 


(1.2) 


This  is  a  tremendous  improvement  over  the  1000  meter  lateral  resolution  of  the  small 
1  meter  physical  radar  antenna.  The  aperture  of  Deff  =  1000  meters  is  called  a 
synthetic  aperture  [Ref.  3]. 

When  SAR  data  are  collected  by  the  radar  system,  the  data  appear  to  be  quite 
unfocused.  In  fact,  the  received  data  look  much  like  random  noise.  As  in  a  hologram, 
the  essential  information  lies  in  the  phase  of  the  received  data,  and  phase-sensitive 
processing  is  needed  to  obtain  a  focused  image  [Ref.  2],  In  the  1950s  and  1960s,  the 
science  of  remote  sensing  was  developing  in  the  civilian  community.  From  the  origins 
of  aerial  photography,  digital  scanners  using  several  optical  frequency  bands  were 
installed  in  aircraft  and  satellites,  and  people  began  developing  uses  for  the  detailed 
wide-area  images  of  the  Earth's  surface  that  were  acquired.  Military  SAR  technology 
was  released  to  the  civilian  community  in  the  1970s,  and  remote  sensing  scientists 
found  that  SAR  images  provided  a  complementary  and  useful  addition  to  their  optical 
sensors  [Ref.  2].  Since  synthetic  aperture  imaging  was  originally  introduced  for  radar 
systems,  this  technology  has  been  used  for  different  military  applications  shown  in 
table  A. 
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•  Automatic  Aircraft  Landing 

•  Mine  detection 

•  Concealed  target  detection  in  foliage 

•Interior  imaging  of  Buildings  in  a  rescue  or  hostage  crisis 

•  Missile  detection  and  tracking 

•  Ground  moving  target  detection  and  tracking 
•Airborne  Reconnaissance 


Table  I.  Military  applications  of  SAR  technology 
In  this  thesis  we  will  focus  on  -  Airborne  Reconnaissance. 

There  are  two  classes  of  SAR  imaging  systems:  spotlight  SAR ;  and  stripmap  SAR.  In 
spotlight  SAR  the  moving  radar  platform  stares  at  a  specific  location  (usually  on  the 
ground)  so  that  at  each  point  in  the  flight  path  the  same  target  is  illuminted  from 
different  directions.  The  locus  of  equal  range  points  is  determined  by  the  intersection 
of  the  radiating  spherical  wave  and  the  ground.  When  the  ground  is  assumed  to  be  a 
plane,  these  points  will  form  a  circle.  If  the  radar  antenna  beamwidth  is  narrow  and 
the  spot  is  sufficiently  far  away,  then  the  data  are  collected  on  small  arcs  that  can 
be  approximated  as  lines.  Consequently,  the  imaging  method  is  mathematically  the 
same  as  that  used  in  Inverse  Synthetic  Aperture  Radar  (ISAR)  (see  Figure  1) 


Figure  1.  In  spotlight  SAR  the  radar  is  trained  on  a  particular  location  as  the  radar 
platform  moves.  In  this  figure  the  equi-range  circles  (dotted  lines)  are  formed  from 
the  intersection  of  the  radiated  spherical  wavefront  and  the  surface  of  a  (flat)  Earth. 
This  Figure  is  taken  from  [Ref.  11] 
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In  contrast,  stripmap  SAR  sweeps  the  radar  beam  along  with  the  platform 
without  staring  at  a  particular  location  on  the  ground  (see  Figure  2).  The  equi-range 
points  are  still  circles,  but  the  data  model  based  on  simple  rotation  is  no  longer 
valid.  In  fact,  since  radar  doesn’t  stare  at  the  same  location,  effective  stripmap 
reconstruction  will  require  wide  radar  beamwith  [Ref.  11]. 


Figure  2.  Stripmap  SAR  acquires  data  without  staring. The  radar  typically  has  fixed 
orientation  with  respect  to  the  platform  and  the  data  are  acquired  as  the  beam 
footprint  sweep  over  the  ground.  This  Figure  is  taken  from  [Ref.  11] 

In  this  thesis  we  will  address  the  construction  of  ISAR  images  using  two  dif¬ 
ferent  approaches:  Fourier  and  Wavelet  transforms.  This  thesis  will  be  organized  as 
following:  In  chapter  two  we  will  derive  the  weak-scatter  far-held  model.  In  chapter 
three  we  will  discuss  the  Radon  transform.  In  chapter  four  we  will  address  the  differ¬ 
ence  between  the  Fourier  and  the  wavelet  transforms.  In  chapter  five  we  will  discuss 
the  algorithms  and  how  to  construct  ISAR  images  using  the  two  approaches.  Finally, 
in  the  sixth  chapter  we  will  summarize  our  conclusions. 
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II. 


THE  RADAR  SCATTERING  MODEL 


A.  WEAK-SCATTER  FAR-FIELD  MODEL 

As  always  we  need  to  start  from  a  mathematical  model,  to  describe  what  is 
happening  in  reality.  In  our  case  we  are  interested  in  the  backscatter  held  from  the 
target  after  being  illuminated  by  an  electromagnetic  wave. 


Figure  3.  Geometry  of  the  radar  imaging  of  a  target.  This  picture  was  taken  from 
[Ref.  11] 

Let’s  assume  that  our  incident  held  is  of  the  form, 

Hlllc  =  H0ei(-kR'r-wt'1  (II.l) 

where  H  =  jt  is  the  constitutive  relation,  B  is  the  Magnetic  Field,  /r  is  the  Magnetic 
permeability  and  k  is  the  wave  number.  If  we  assume  the  target  to  be  a  perfect 
conductor  (current  just  on  the  surface)  and  that  its  surface  is  smooth,  then  after 
some  algebra  we  can  define  the  scatter  held  from  a  target  as  follows,  [Ref.  11] 

ikHr.pi(2kR-wt'>  r  ~  , 

Hscatt(R,  t:k)=  ^  J '  PkMx/)ei2kRx  dx'dy'dz'  (II.2) 


This  is  called  the  weak-scatterer  far-held  model.  The  function  pk^xi\  is  dehned  as 
a  generalized  scatterer  density  function  which  will  generally  depend  on  k  and  R.  It 
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can  be  shown  that  the  integrand  of  equation  (II. 2)  is  the  far  held  pattern  F(9,<j)). 
This  is  why  (II. 2)  is  called  the  weak-scatter  far-held  model.  We  are  interested  in  this 
expression  because  from  Hscatt  we  can  extract  information  about  the  target  of  interest 
by  estimating  various  functions  closely  related  to  pk  p.  For  our  purposes,  (the  imaging 
case),  it’s  useful  to  dehne  a  cartesian  coordinate  system  hxed  to  the  target.  This 
coordinate  system  will  move  accordingly  with  the  target  as  it  maneuvers.  Another 
assumption  is  that  the  target  moves  as  a  rigid  body  and  that  its  instantaneous  axis 
of  rotation  has  a  component  k  _L  R.  To  describe  rotations  of  the  target  ,  we  need 
to  introduce  a  new  reference  frame  by  defining  the  directions  i  and  j  in  terms  of  the 
rotation  angle  6  by  i  ■  R  =  sinO  and  9  by  j  •  R  =  cos6 :  these  are  the  cross-range  and 
down-range  directions,  respectively  (see  Figure  3). 

Now  using  the  new  coordinates,  equation  (1 1. 2)  can  be  written  as 


H  scatt^-F^i  t  ^  ) 


ikH0ei(-2kR~wV 

(2tt)3R 


Pk, R(X'Y ei2k(y'cos9-x'sine)  dx'dy'dz' 


(II.3) 


and  x'  =  x'i  +  y'j  +  z'k.  This  is  the  standard  target  model  used  in  radar  imaging 
analysis.  We  know  that  pk  p  depends  on  k  and  R ,  this  dependence  which  is  often 
treated  as  insignificant  over  the  span  of  R,  and  the  ranges  of  k  obtainable  from  con¬ 
ventional  radar  systems.  This  dependence  is  typically  ignored.  This  approximation 
allows  us  to  exploit  equation  II. 3  as  a  Fourier  transform  relationship  between  H scatt 
and  p  =  pkp.  The  Founder  transform  is  very  convenient  for  practical  object  function 
estimation,  but  this  approximation  turns  out  to  cause  significant  problems  in  image 
interpretation  [Ref.  11]. 

Radar  targets,  especially  man  made  targets,  can  be  considered  as  a  collection 
of  point  scatterers.  The  behavior  of  these  scatterers  has  a  vast  variety  of  reflecting 
and  back-scattering  types.  They  can  be  edges,  corners,  circles,  cavities,  squares  (see 
Figure  5) 
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CORNER 

DIFFRACTION 


Figure  4.  Scatter  mechanism. 


This  picture  was  taken  from  [Ref.  12] 


Figure  5.  Different  back-scattering  behavior 


B.  SCATTERING  FROM  A  SPHERE 

To  better  understand  how  scattering  depends  on  frequency  let’s  introduce  the 
simple  example  of  scattering  from  a  sphere.  Assuming  a  plane  harmonic  wave  of 
frequency  u>  (with  k  =  and  traveling  along  the  z  axis  is  incident  on  a  perfectly 
conducting  sphere  of  radius  a.  Let’s  determine  the  physical  optics  contribution  to  the 
scattered  held  and  the  associated  scattering  cross-section  of  the  sphere  with  radius 
r=l. 


ikH0ei{-2kR~wt ] 
2v tR 


'  Rn<o 


R  ■  hei2kR'x'  dS' 


(II.4) 
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where  f  =  h  =  r  sin  9  cos  4>i  +  r  sin  9  sin  <j)j  +  r  cos  9k  (in  spherical  coordinates) 


R  =  k 

(II.5) 

h  ■  R  —  r  cos  9  =  cos  9 

(II.6) 

R-  x'  =  R  -  ah  =  a  cos  9 

(II.7) 

since  just  half  of  the  sphere  is  visible 


R  he1 


dS'  =  2iraz  I  '  cos9ei2kacosG  sin  9d9 
o 


J  Rh<  o 

Making  now  the  substitution  cos  9  =  u  and  —  sin  9d9  =  du  we  obtain 


(II.8) 


R  he 1 


dS’  -  2iraz  I  '  uei2kaudu 
o 


JRh<  o 

After  some  algebra  we  find  the  result  of  this  integral  to  be 


(IL9) 


ttR.*' 9 1  —  ei2fc“  +  2emaika 


iRh<  o 


R  •  he dS '  =  2 7raz 


i24k2a2 


therefore 


(11.10) 


Hsra / /  ( R‘- 1 .  &) 

which  reduces  to 


ikH0ei(2kR~wt)  o  1  -  ei2fca  +  2 ei2kaika 
-2ira'J  — 


2vr  R 


4  k2a2 


(11.11) 


•  rj  _  pi(2kR—wt) 

Hscatt(R,  t:k)=  °  A1  - - (1  -  ei2fca  +  2 ei2kaika) 


4  kaR 


We  define  the  radar  cross-section  a  by 


(11.12) 


a  \Hsrntt\2  1 


7 ra2  I H o  1 2  7ra2 


where  H 0  is  the  incident  field  [Ref.  11]. 


(11.13) 
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Figure  6  shows  the  basic  regions  characterizing  the  radar  cross  section  of  a 
sphere.  The  optical  region  (often  called  the  Fraunhofer  region )  applies  when  >  10. 
In  this  region,  the  cross-section  is  independent  of  frequency.  At  point  A  of  Figure  6 
we  have  a  perturbation  ,  the  largest  one,  that  would  be  4  times  higher  than  the  cross 
section  found  in  the  optical  region.  A  minimum  also  occurs  at  point  B  where  the 
cross  section  would  be  0.26  times  the  value  in  the  optical  region  formula.  This  area 
is  known  as  Mie  or  resonance  region.  There  is  another  region  where  the  diameter 
of  the  sphere  is  smaller  than  the  wavelength:  the  Rayleigh  region.  As  we  can  see, 
the  cross-sectional  is  smaller  than  in  the  other  regions  and  this  happens  because  the 
wavelength  is  larger  than  the  sphere,  therefore,  the  scattered  held  in  the  direction  of 
the  sensor  is  very  weak.  Rayleigh  scattering  is  important  in  weather  radar. 

In  our  case  we  will  be  in  the  optical  region  because  the  bandwidth  to  be  used  in 
the  test  cases  is  from  8.75GHz-9.25Ghz.  The  wavelength  is  of  the  order  of  the  order 
of  3cm,  so,  sufficient  to  detect  parts  of  airplanes,  ships,  cars  .  As  we  saw  in  the 


Figure  6.  Coss-section  plot  on  a  log-log  scale  as  a  function  of  kr 
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weak-scatter  far-field  model,  because  we  are  in  the  optical  region  we  will  have  some 
artifacts  in  the  images  crated  by  corner  reflections,  inlet  engines  and  other  types  of 
undesired  reflections  from  other  geometrical  shapes  of  the  target. 


C.  IMAGING  FROM  WEAK-SCATTERER  FAR-FIELD 
DATA  USING  FOURIER  ANALYSIS 

We  can  see  from  the  weak-scatter  far-held  model  of  equation  11.14  that  the 
data  is  similar  to  a  Fourier  transform.  In  fact  it  is. 


{FPe}{u)  = 


1 


r»oo  POO 


p(x.  y)e<“{x'in‘-yc°,,)ixdy 


(11.14) 


-oo  J  —oo 


In  this  approach  we  will  take  data  over  a  small  aperture  A 9  with  0  <  9  <  A9  <C  1. 
The  small  angle  approximation  sinO  ~  9  and  cos9  ~  1  implies 

1 


r»oo  poo 


{ FPe }  ~  M- 


p{x,y)elu(yXd  yS> dxdy 


(11.15) 


-oo  J  — OO 


a/27 t  J-c 

As  we  can  see  now  the  exponential  in  the  equation  11.15  has  linear  argument,  since 
kx  =  uj9  and  ky  —  u  and  so  this  equation  can  now  be  inverted  using  Fourier  analysis: 
From  the  frequency  domain  data  set{Fpe}  i— >  P(kx,  ky),  can  be  recovered  as  p(x,y) 


i  POO  POO 

P{x.y)=^—  /  P(kx,ky)e*k*x-k^dkxdky  (11.16) 

For  reference  purposes  we  switch  x  and  y  in  the  exponential  argument  from  the 
weak-scatter  far-held  model  and  equation  11.15  we  label  the  cross  range  coordinate 
as  x  and  the  down-range  coordinate  as  y. 

The  data  set  mapping  {Fpg}  >  P(kx,ky )  is  based  on  the  transformation  kx  =  uj6 
and  ky  =  lu.  If  we  consider  ky  =  uj  to  be  a  radial  direction  in  the  data  space,  then 
it’s  easy  to  see  that  kx  =  lo9  =  Ky9  is  an  arclength  for  hxed  ky.  As  was  mentioned 
before,  we  are  assuming  small  angles  because  our  data  set  will  be  collected  over  a 
small  angle  grid  in  polar  coordinates  (see  Figure  7). 

We  know  that  equation  (11.16),  on  the  other  hand,  is  a  two-dimensional  Founder 
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transform  formula  appropriate  to  data  defined  on  a  rectangular  grid  —  that  is,  in 
this  approach  the  data  is  processed  on  a  rectangular  grid  but  collected  on  a  polar  one 
[Ref.  11]. 


Figure  7.  Data  space  used  in  the  image  recovery  equation  (11.16).  Here  ky  —  u  is  a 
polar  direction  and  ky  =  u6  is  an  arclength  (for  each  u).  This  picture  was  taken  from 
[Ref.  11] 


Figure  8.  Radar  Signals  capture  method.  This  picture  was  taken  from  [Ref.  12] 
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When  uj  is  large  and  the  aperture  is  sufficiently  small,  the  polar  grid  is  pretty 
close  to  rectangular  in  shape.  But  in  the  general  case,  the  measured  data  must  be 
interpolated  to  a  rectangular  grid  before  equation  11.16  can  be  applied  (an  illustration 
is  in  Figure  9).  In  the  radar  community,  this  additional  preprocessing  step  is  known 
as  focussing  [Ref.  11]. 


polar  data  grid 


rectangular  / 
interpolating  grid 


o  —  measured  data 
•  —  interpolated  data 


Figure  9.  Data  interpolation  from  a  polar  grid  to  a  rectangular  grid.  This  picture 
was  taken  from  [Ref.  11] 

One  of  the  relevant  aspects  is  the  resolution  of  the  image:  equation  (11.16)  gives 
us  the  means  to  evaluate  the  resolution  of  the  reconstructed  image.  This  equation  is 
the  Fourier  domain  version  of  a  convolution  with  an  imaging  kernel  determine  by  the 
inverse  Fourier  transform  of  k(kx,  ky)  -  e^kx+ky\ 

Since  we  deal  with  finite  bandwidth  signals,  the  spatial-domain  version  of  this  imaging 
Kernel  is  a  product  of  sine  functions  and  the  width  of  the  central  lobes  provides  an 
estimate  of  the  resolution.  The  dimensions  Ax  and  Ay  of  a  resolution  cell  can  be 
estimated  as 


2n  2i t 


(11.17) 


A  kx  uA  9 


2tt  2tt 


(11.18) 
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where  uJ  is  the  average  frequency,  A 6  represents  the  aperture  size,  and  A uj  denotes 
the  bandwidth  [Ref.  11].  We  can  see  that  the  down-range  resolution  depends  on  the 
bandwidth  and  the  cross-range  resolution  depends  on  the  aperture  size,  scaled  by 
frequency.  Don’t  forget  we  are  using  the  small  angle  approximation  and  it’s  because 
of  that  assumption  that  this  analysis  is  valid.  This  method  will  be  used  on  test  cases 
in  the  following  chapters  and  analyzed  concerning  the  pros  and  cons  of  this  technique. 


Figure  10.  Standard  Fourier  image.  This  Figure  was  taken  from  [Ref.  13] 
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III. 


RADON  TRANSFORM 


A.  PRELIMINARY  DISCUSSION 

The  objective  of  our  work  is  to  obtain  range  and  cross  range  target  information 
from  High  Range  Resolution  (HRR)  data.  To  understand  how  to  achieve  this  goal, 
let’s  examine  a  data  model  for  an  impulsive  signal  (ideal  short  pulse  in  time  domain) 

/OO  /*oo 

/  p(t',  u')S(t  —  r/)eu'(r-T  ^ dr'du '  =  p(r)  (III.  1) 

-oo  J  — OO 

The  HRR  data  are  formed  as  the  convolution  between  the  first  range  argument  of 
p(r,  v)  with  the  measurement  kernel  (idealized  by  the  delta  function).  In  our  case  we 
want  to  distinguish  range  profiles  by  the  target  aspect  angle  at  which  the  data  are 
collected.  Therefore, 


/OO  /»oo 

/  Pe(j',  u')5(t  -  Tl)e%v^r~"r'^ dr' dv'  =  pe(r)  (III. 2) 

-oo  J  —  oo 

since  the  scattering  density  function  p(r' ,  o')  will  depend  on  target  orientation  (see 
Figure  11  ). 


radar 


Figure  11.  Geometry  of  an  HRR  profile  at  an  angle  6.  This  picture  was  taken  from 
[Ref.  11] 
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But  the  function  pe(r',  v')  is  just  a  rotated  version  of  the  stationary  target.  If 
we  assume  the  coordinate  system  x  —  y  is  fixed  to  the  target  and  set  the  scattering 
density  function  in  this  coordinate  system  by 

p{x,y)  =  pe=0{r',v') 

then  for  a  general  6  it’s  easy  to  see  that, 

Pq{t u')  =  p(—r'sind  +  u'cosd ,  t'cosO  +  v'sind) 


Using  this  result  in  equation  III. 3  yields 
Ve(r,u)  = 


r*00  POO 


'  — oo  J  — oo 


p{—r'sinQ  +  u'cosd,  t'cosO  +  u'sind)5{r  —  r')etl'(-T  dr'du '  =  p(r) 

(HI-3) 

Under  the  change  in  variables 


t  =  —xsind  +  ycosQ 


v'  =  +xcos8  +  ysinO 


equation  III. 3  becomes 


f »oo  /»oo 


Ve(r,v)  = 


'  — oo  J  — oo 

POO  POO 


p(x,  y)S(r  +  xsind  -  ycosd)eiv{T+xsin(,-ycose) dxdy 

p{x ,  y)6(r  +  xsind  —  ycos6)dxdy  (III. 4) 


'  — oo  «/  —  OO 


As  we  can  see,  this  is  just  a  line  integral  of  p(x,y)  along  the  line  L(r\d)  defined 
by  r  +  xsind  —  ycosd  —  0  (because  of  the  (5-function  in  the  integrand).  This  line  is 
illustrated  in  Figure  11.  Note  also  that  the  argument  of  the  complex  exponential  is 
zero  and  so  doesn’t  contribute  to  the  integral. 

The  r}e(T,n)  in  equation  III. 4  is  independent  of  v  and  the  double  integral 
restricted  to  the  line  L(r;  d)  or  (the  equivalent  line  integral  itself)  is  known  as  the 
Radon  Transform  of  p(x,y)[ Ref.  11], 


{Rp}(r,d)  = 


POO  POO 


p(x,  y)6(r  +  xsind  —  ycosd)dxdy 


(III.5) 


'  — oo  J  — oo 
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The  factor  S(t  +  xsinO  — ycosO )  is  zero  everywhere  except  where  its  argument  is  zero, 
which  is  along  the  straight  line  xsinO  —  ycosd  =  r.  The  straight  line  L  represents  a 
slit  when  it  is  at  a  perpendicular  distance  r  from  the  origin  and  inclined  at  an  angle 
9  to  the  y-axis.  If  9  is  kept  fixed,  say  a  value  6\  while  r  is  varied,  then  the  integral 
(i?p}(r,  0\)  constitutes  the  projection  of  the  density  function  p(x,y )  onto  the  line 
9  =  #i  as  a  function  of  r. 

A  practical  computational  method  for  inversion  can  be  arrived  at  by  Fourier  trans¬ 
forming  the  Radon  integral  equation,  finding  a  method  of  solution,  and  then  trans¬ 
forming  the  steps  to  end  up  with  a  data-plane  operations  in  which  numerical  Fourier 
transformation  is  actually  dispensed  with.  The  inversion  procedure  derives  from  a 
remarkable  relation  that  exists  between  Fourier,  Abel,  and  Hankel  transforms  and 
from  a  generalization  known  as  the  Projection-Slice  Theorem[Ref.  5].  This  theorem 
is  discussed  in  the  next  section. 

B.  BACKPROJECTION 

For  every  9  define  |i?p}(r,  9)  =  ~pe{r)  of  equation  11  to  be  a  function  of  one 
variable  r  then  r  for  each  9.  The  Fourier  transform  of  pg(r)  can  be  computed  as 


Fpe{u) 


\/2 W_ 


J  —  OO 

r*oo 


'  — OO 

poo 


Pe(r)eiuTdr 

p(x,y)S(r  +  xsin9  —  ycos9)elbJT  drdxdy 
p(x,y)ei^xsine-ycose)dxdy 
p(x,y)ei{k*x-kyy)dxdy 


(III. 6) 


assuming  kx  =  —ujsin9  and  ky  =  ujcos9. 

We  observe  that  equation  III. 6  corresponds  to  a  2D  Fourier  transform  of  the 
density  function  p(x,y).  As  is  well  known,  the  2D  Founder  transform  is  given  by 


1  poo  poo 

P(kx,  ky)  =  —  /  p(x,  y)ei^x-k^dxdy 

2vr 


(III.7) 
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and  so  we  can  conclude 


P(—usinQ,ujcos6 )  =  F  pe  (o>) 

v27T 


(III.S) 


Equation  III. 8  is  known  as  the  Projection- Slice  Theorem ,  and  allows  us  to  determine 
p(x,  y)  from  the  inverse  Fourier  transform  of  F{kx ,  ky).  We  can  use  this  result  to  build 
an  inverse  Radon  transform  for  p(r,  6).  Using  the  change  of  variables  (. kx ,  ky)  i— >  (a;,  6) 
where  kx  =  —usind  and  ky  =  locosO ,  we  obtain 


P(X,V)  =  ^  fr  P(kx,ky)^-^dkxdky 


1 


r»27T  POO 


)  d(kx,  ky 


/  /  Z7—  ,\„—iuj(xsinO—ycos6)'^VVx'>'y'yJj, 

Jo  L  Fp,{")e 

p2tt  poo 


l 


(27t)3/2  Jq 

j  p2tt  poo 


Fpe  (^e-Mxsine-ycosO)  |  w  | 


u\ Fpe{u)  )  eM-xsind+ycose)dLud6  (III. 9) 


(2tt)3/2  J0  J 

The  inner  integral  of  equation  III. 9  is  just  the  inverse  Fourier  transform  of  the 


product 


Q(uj)  =  w  x  {Fpe(cj)} 
The  inverse  Fourier  transform  of  Q( cu)  is 


1  f°° 

q(t )  =  -f=  /  Q^)ewt 

V  J — oo 

and  so  the  equation  III. 9  can  be  written  as 


r*27T 


p(x,y)  =  —=  /  q(—xsin6  +  ycos6)d0 
v2vr  Jo 


where  after  substitution  of  Fpe(uj),  we  obtain 


1  POO  POO 

^=2~J  I 


'  — oo  */  — oo 


(III.10) 


(III.ll) 
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Equation  III. 10  is  called  the  the  backprojection  operator  and  the  filtered  backprojection 
algorithm  combines  equations  III.  10  and  III.  11  to  form  the  an  inverse  Radon  transform 
P(T,6 0  =  p(r)  ->  p(x,y). 


a.  Using  3  views 


b.  Using  many  views 


Figure  12.  Backprojection  reconstructs  an  image  by  taking  each  view  and  smearing 
it  along  the  path  it  was  originally  acquired. The  result  image  is  a  blurry  version  of  the 
correct  image.  This  picture  was  taken  from  [Ref.  15] 


filtered 
view  3 


filtered  view  1 

i  1 

zs:x . 
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a.  Using  3  views 


b.  Using  many  views 


Figure  13.  Filtered  backprojection  reconstructs  an  image  by  filtering  each  view  before 
backprojection.  This  removes  the  blurring  seen  in  simple  backprojection,  and  results 
in  a  mathematically  exact  reconstruction  of  the  image.  Filtered  backprojection  is  the 
most  commonly  used  algorithm  for  computed  tomography  systems.  This  picture  was 
taken  from  [Ref.  15] 
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IV.  TIME-FREQUENCY  METHODS 


A.  ISAR  IMAGING  USING  TIME  FREQUENCY  TECH¬ 
NIQUES 

The  approach  to  Inverse  synthetic  aperture  radar  (ISAR)  is  exactly  the  same 
as  for  SAR,  but  in  ISAR  we  exploit  the  target’s  motion.  The  goal  is  to  generate 
images  of  moving  targets  in  the  range-Doppler  plane.  Usually  the  Fourier  transform 
is  used  to  compute  the  image,  but  because  of  the  unpredictable  movement  of  the  tar¬ 
get  the  Doppler-frequency  can  display  nonuniform  time  variations.  To  overcome  this 
time- varying  Doppler  spectrum  we  investigate  the  use  of  a  time-frequency  transform. 

In  this  thesis  we  will  examine  the  performance  of  two  different  transforms,  the 
short  time  Fourier  transform  itself  and  the  wavelet  transform  to  obtain  ISAR  images. 
Time-frequency  techniques  are  used  to  overcome  limitations  of  the  standard  Fourier 
transform.  For  example,  during  the  collection  of  data  the  scattering  points  must  re¬ 
main  in  the  same  range  cell:  if  this  doesn’t  happen,  the  image  will  be  blurred.  This 
occurs  when  scatterers  move  out  of  their  range  cells  making  our  Doppler  information 
non-linear.  Since  the  Fourier  transform  requires  linear  sampling,  we  can’t  compensate 
for  this  effect.  To  obtain  a  focussed  radar  image  using  the  Fourier  transform,  motion 
compensation  algorithms  are  necessary  to  adjust  for  scatterer  migration.  But  this 
is  often  not  enough  because  the  spectrum  can  still  display  non-linear  behavior.  So, 
Doppler  tracking  and  other  techniques  must  be  applied.  When  a  target  is  moving 
smoothly,  these  standard  techniques  are  normally  good  enough  to  focus  our  image  and 
we  can  still  use  the  Fourier  transform.  However,  when  the  target  exhibits  unexpected 
motion,  such  as  rapid  maneuvers  and  acceleration  about  an  axis,  the  standard  motion 
compensation  algorithms  are  not  enough  to  focus  the  image.  In  this  situation,  more 
advance  algorithms  for  motion  compensation  are  needed  [Ref.  16].  However,  such 
advanced  algorithms  are  not  always  effective  against  the  image  blurring,  because  the 
scatterers  may  still  drift  out  of  their  range  cells  and,  as  a  consequence,  their  Doppler 
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frequency  shifts  become  nonlinear  in  time.  Therefore  the  image  will  still  come  out 
blurred  after  the  Fourier  transform. 

As  a  matter  of  fact,  the  Fourier  transform  limitations  can  be  mitigated  using 
the  Short  Time  Fourier  Transform  (STFT).  As  we  explained  before,  the  Doppler 
information  is  nonlinear  with  time,  therefore  a  good  method  to  solve  the  smeared 
Fourier  frequency  spectrum  is  to  apply  high-resolution  time-frequency  transform  to 
the  processing.  Using  this  method,  we  no  long  need  the  motion  compensation  algo¬ 
rithms  because  we  are  no  longer  using  the  Fourier  transform.  This  way,  unpredictable 
Doppler  frequency  shifts  can  be  overcome  without  using  sophisticated  motion  com¬ 
pensation  algorithms.  The  use  of  the  STFT  will  be  the  first  to  be  applied  to  our  case. 
In  the  same  way,  another  frequency  analysis  tool  can  be  applied  using  the  Continuous 
Wavelet  transform  (CWT).  The  wavelet  approach  offers  an  alternative  path  for  sig¬ 
nal  processing.  Essentially,  the  CWT  is  the  same  as  the  STFT  except  that  any  basis 
functions  can  be  used  (only  harmonic  functions  of  constant  amplitude  and  phase  are 
used  by  the  STFT  [Ref.  17].  Another  difference  between  the  CWT  and  the  STFT 
is  the  window  size:  the  STFT  uses  one  window  analysis  with  the  same  width,  while 
the  CWT  uses  different  windows  for  different  frequencies.  (Short  windows  for  high 
frequencies  and  long  windows  for  low  frequencies.) 

1.  Short  Time  Fourier  Transform  vs.  Continuous 

Wavelet  transform 

In  the  standard  FT  we  assume  signals  of  infinite  length  with  stationary  spec¬ 
trum,  thus  making  it  not  well  suited  to  the  detection  of  time  varying  spectra.  The 
STFT  on  the  other  hand  computes  the  FT  of  the  signal  on  a  window  sliding  in  time, 
thus  capturing  the  non-stationary  behavior.  The  STFT  is  defined  in  the  following 
way 

/OO 

s(t)w(t  -  t)e~iUTdr  (IV.l) 

•OO 

The  window  w(r  —  t)  is  centered  on  the  time  variable  t  and  then  shifted  the  window 
through  out  the  signal  to  obtain  the  information  about  the  frequency  behavior  in  time 
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(see  Figure  14)  .  As  far  as  its  application  to  imaging  is  concerned  STFT  is  affected 
by  the  uncertainty  principle  by  which  the  signal  cannot  be  resolved  in  both  time  and 
frequency  (see  Figure  15  and  Figure  16  ).  Also,  the  STFT  is  only  a  window-averaged 
frequency  spectrum  and  not  an  exact  estimation  of  the  instantaneous  frequency  [Ref. 
18]- 


Figure  14.  Short  time  Fourier  Transform.  This  figure  was  taken  from  [Ref.  10] 


w(t)  \W(F)\ 


Figure  15.  Time/frequency  uncertainty  of  a  rectangular  window.  This  figure  was 
taken  from  [Ref.  10] 
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Figure  16.  Basis  functions  and  time/frequency  resolution  for  STFT.  This  figure  was 
taken  from  [Ref.  10] 


The  best  compromise  is  to  have  different  resolutions  at  different  frequencies 
by  adapting  the  window  length  in  time.  This  is  what  the  CWT  tries  to  accomplish. 


CWT(t,  a)  =  /  x(t)a°'5ip(a(t  —  r))dt  (IV.2) 

J  — OO 

As  we  saw  before,  the  STFT  has  a  fixed-length  window  which  is  responsible 
for  the  principal  limitation  of  time  and  frequency  resolution  due  to  the  uncertainty 
principle.  Either  we  resolve  in  time  or  we  resolve  in  frequency.  In  order  to  circumvent 
this  issue,  we  will  observe  the  signal  using  different  window-lenghts  because  the  CWT 
allows  it.  Normally  we  are  interested  in  the  higher  frequencies  due  to  rapid  function 
variation  as  a  function  of  time.  Evidently,  it  makes  some  sense  to  use  a  long  window 
for  low  frequencies  and  a  short  window  for  high  frequencies  to  get  high  resolution 
over  all  spectrum  specially  at  higher  frequencies. 
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Figure  17.  Basis  functions  and  time/frequency  resolution  for  CWT.  This  figure  was 
taken  from  [Ref.  10] 


Figure  18.  Basis  functions  and  time/frequency  resolution  of  the  STFT  and  the  CWT. 
The  tiles  represent  for  a  given  basis  functions  its  concentration  in  the  time-frequency 
plane,  a)  Coverage  of  the  time-frequency  plane  for  the  STFT.  b)  for  the  CWT  c) 
Corresponding  basis  functions  for  the  STFT.  d)  for  the  CWT.  This  figure  was  taken 
from  [Ref.  19] 
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a) 


Constant  Bandwidth  (STFT  Case) 


Figure  19.  Division  of  the  frequency  domain  a)  for  the  STFT  (uniform  coverage)  and 
b)  for  CWT  (logarithmic  Coverage).  This  figure  was  taken  from  [Ref.  19] 


B.  EXAMPLE  OF  A  CHIRP  USING  STFT  VS.  CWT 

We  present  in  Figures  20  and  21  a  spectrogram  and  a  scalogram  from  a  chirp 
signal  of  length  2000  samples,  where  the  frequency  grows  continuously  and  linearly. 
The  spectrogram  is  obtained  using  a  Gaussian-like  window  and  the  scalogram  is 
obtained  with  a  Morlet  wavelet  (a  complex  sinusoid  windowed  with  a  Gaussian  en¬ 
velope).  In  this  way,  we  can  see  that  for  the  STFT  the  width  of  the  frequency  line  is 
always  the  same  as  a  result  of  the  constant  window  width  and  the  uniform  coverage 
in  the  frequency  domain.  However,  the  scalogram  with  nonlinear  (logarithmic)  scale 
results  in  a  curve  instead  of  a  line.  Also,  we  can  see  that  the  width  of  the  line  is 
thiner  for  high  frequencies,  a  result  that  reflects  the  long  duration  windows  in  time 
for  higher  frequencies. 
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Figure  20.  Chirp’s  Spectrogram 


Figure  21.  Chirp’s  Scalogram 
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V. 


ISAR  IMAGING  USING  STFT  AND  CWT 


A.  BACKHOE  DATA  DOME 

To  further  examine  the  potential  of  using  time-frequency  methods  to  produce 
ISAR  images,  we  examined  simulated  data  based  on  the  Air  Force  Research  Labora¬ 
tory  (AFRL),  Backoe  Data  Dome,  Version  1.0. 

1.  Backoe  Data  Dome  Version  1.0 

The  Backhoe  Data  Dome,  Version  1.0  consists  of  simulated  wideband  (7-13 
GHz),  full  polarization,  complex  backscatter  data  from  a  backhoe  loader  generated 
using  a  Computer-Aided-Design  (CAD)  model  in  free  space  as  shown  in  Figure  22. 
The  MATLAB  program  plotfacet.m  available  from  the  Air  Force  Research  Laboratory 
Sensor  Data  Management  System  was  used  to  plot  all  3-D  graphics  of  the  backhoe 
loader. 


Figure  22.  CAD  Model  of  Backhoe  Loader 


The  backscatter  data  consists  of  over  1,125  Matlab  binary  files  (*.mat  files) 
in  which  each  file  has  the  k  -  space  data  from  with  elevation  between  18  and  42 
degrees  and  Azimuth  between  66  and  114  degrees.  The  filenames  are  indicative  of 
the  data  they  contain.  backhoe-el042-azl09.mat,  for  example,  contains  the  data  with 
elevation  between  29.5  and  30.5  and  azimuth  between  89.5  and  90.5.  The  viewing 
hemisphere  consists  of  14  samples  per  degree  in  both  azimuth  and  elevation.  There 
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is  one  sample  every  11.75  MHz  over  the  6  GHz  bandwidth  and  full  polarization 
information;  i.e.  horizontal  linear  polarization  on  transmit  and  receive  (HH),  vertical 
linear  polarization  on  transmit  and  receive  (VV),  and  horizontal  linear  polarization 
on  transmit  and  vertical  linear  polarization  on  receive  (HV). 


Figure  23.  Backhoe  Data  Dome  Representation  in  k-Space 

With  the  methods  described  in  the  preceding  Chapters,  the  STFT  and  CWT 
will  be  applied  to  the  3-D  k-space  data  for  the  different  test  cases,  with  the  following 
setup  in  the  following  orientations  and  using  the  specified  polarization: 

Test  Case  1:  Data  collected  over  the  following  interval  (horizontal  Path): 
Elevation=42deg  and  95deg  <  Azimuth  <  102deg. 

Test  Case  2:  Data  collected  over  the  following  interval  (curve  Path):  38<  El¬ 
evation  <  42deg  and  95deg  <  Azimuth  <  99deg. 
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In  all  cases  the  center  frequency  chosen  was  9.25GHz  with  a  Bandwidth  of 
500  MHz.  We  chose  500  MHz  because  it’s  a  more  realistic  case,  and  is  representative 
of  existing  radars  in  the  world.  With  respect  to  Polarization,  this  will  not  be  very 
important  in  our  simulation  but,  of  course,  Polarization  matters  in  the  case  we  have 
a  Radar  with  cross-polarization  capabilities.  In  this  case,  the  amount  of  information 
available  to  reconstruct  the  image  will  be  greater  than  with  just  vertical  or  horizontal 
polarization.  Cross-pol  images  will  be  briefly  analyzed  below. 

B.  ALGORITHM  FOR  ISAR  IMAGING  USING  TIME- 
FREQUENCY  METHODS 

After  the  complex  data  is  extracted  from  our  database,  we  construct  a  matrix 
in  which  the  N  columns  have  the  range  information  and  the  M  rows  have  the  cross¬ 
range  information. 

step  1.  We  do  the  Fast  Fourier  Transform  of  the  rows  to  obtain  our  range  pro¬ 
files.  For  each  row  X{1)  =  FFT[xi(n)\,n  =  1,2,...,  A;  i  =  1,2, 

step  2.  Compute  the  STFT  or  CWT  using  different  windows  type,  sizes  and  scales. 
Y(m)  =  STFT[yi(n )]  or  Y(m)  =  CWT[yi(n)],  n=l,2 ,...,]V;  i=l,2,...,M 

step  3.  Find  the  absolute  value  for  each  entry  of  the  matrix  and  plot  it. 

Step  1  to  3  form  the  algorithm  for  ISAR  imaging  using  the  two  time-frequency 
methods:  first  Short  Time  Fourier  Transform ;  and  second  the  Continuous  Wavelet 
Transform. 
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C.  TEST  CASE  1  :  95°  <  6 azimuth  <  102°,  cf) elevation  —  42° 


Cross-range 


Figure  24.  ISAR  image  using  STFT  with  a  hamming  window  with  20  points  of  size 
and  50%  of  overlaping  (Horizontal  Polarization)  -  ( 9th  frame) 


Cross-range 


Figure  25.  ISAR  image  using  STFT  with  a  hamming  window,  size  40  points  and  50 
percent  of  overlapping  (Horizontal  Polarization)  -  ( 2th  frame) 


We  can  already  observe  that  when  we  increase  the  window  size  the  image  is  going  to 
display  better  resolution. 
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Cross-range 


Figure  26.  ISAR  image  using  STFT  with  a  hamming  window,  size  50  points  and  50 
percent  of  overlapping  (Horizontal  Polarization)  -  ( 3th  frame) 


Cross-range 


Figure  27.  ISAR  image  using  STFT  with  a  hamming  window,  size  60  points  and  50 
percent  of  overlapping  (Horizontal  Polarization)  -  ( 2th  frame) 
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Down-range 


Cross-range 


Figure  28.  ISAR  image  using  STFT  with  a  hamming  window,  size  80  points  and  50 
percent  of  overlapping(Horizontal  Polarization)  -  ( 1th  frame) 


Cross-range 


Figure  29.  ISAR  image  using  CWT  with  a  mexican  hat  window,  using  a  scale  of  1:5 
(Horizontal  Polarization) 


Cross-range 


Figure  30.  ISAR  image  using  CWT  with  a  mexican  hat  window,  using  a  scale  of  3:6 
(Horizontal  Polarization) 


Cross-range 


Figure  31.  ISAR  image  using  CWT  with  a  mexican  hat  window,  using  a  scale  of 
10:13  (Horizontal  Polarization) 


Cross-range 


Figure  32.  ISAR  image  using  CWT  with  a  Morlet  window  (Horizontal  Polarization) 
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Cross-range 


Figure  33.  ISAR  image  using  CWT  with  a  Daubechies  2  window  (Horizontal  Polar¬ 
ization) 


D.  TEST  CASE  2  :  38<  ELEVATION  <  42DEG  AND 

95DEG  <  AZIMUTH  <  99DEG. 

When  we  use  STFT  we’ll  observe  the  resolution  increasing  when  we  increase  the 
window’s  size.  When  we  use  CWT  we  can  observe  that  when  we  change  scale  we  can 
loose  information  or  improve  the  resolution  of  the  STFT  approach. 


Figure  34.  ISAR  image  using  STFT  with  a  hamming  window,  size  20  points  and  50 
percent  of  overlapping  (Horizontal  Polarization)  -  ( 8th  frame) 
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Cross-range 


Figure  35.  ISAR  image  using  STFT  with  a  hamming  window,  size  40  points  and  50 
percent  of  overlapping  (Horizontal  Polarization)  -  ( 3th  frame) 


Cross-range 


Figure  36.  ISAR  image  using  STFT  with  a  hamming  window,  size  50  points  and  50 
percent  of  overlapping  (Horizontal  Polarization)  -  ( 2th  frame) 
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Cross-range 


Figure  37.  ISAR  image  using  STFT  with  a  hamming  window,  size  60  points  and  50 
percent  of  overlapping  (Horizontal  Polarization)  -  ( 2th  frame) 


Cross-range 


Figure  38.  ISAR  image  using  STFT  with  a  hamming  window,  size  80  points  and  50 
percent  of  overlapping  (Horizontal  Polarization)  -  ( 1th  frame) 
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Cross-range 


Figure  39.  ISAR  image  using  STFT  with  a  hamming  window,  size  100  points  and  50 
percent  of  overlapping  (Horizontal  Polarization)  -  ( 1th  frame) 


Cross-range 


Figure  40.  ISAR  image  using  CWT  with  a  mexican  hat  window,  using  a  scale  of  1:5 
(Horizontal  Polarization) 
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Cross-range 


Figure  41.  ISAR  image  using  CWT  with  a  mexican  hat  window,  using  a  scale  of  3:6 
(Horizontal  Polarization) 


Cross-range 


Figure  42.  ISAR  image  using  CWT  with  a  mexican  hat  window,  using  a  scale  of 
10:13  (Horizontal  Polarization) 


Cross-range 


Figure  43.  ISAR  image  using  CWT  with  a  Morlet  window  (Horizontal  Polarization) 
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Cross-range 


Figure  44.  ISAR  image  using  CWT  with  a  Daubechies  2  window  (Horizontal  Polar¬ 
ization) 


Cross-range 


Figure  45.  ISAR  image  using  STFT  with  a  hamming  window,  size  40  points  and  50 
percent  of  overlapping  (Cross  Polarization)  -  ( 3th  frame) 
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VI. 


CONCLUSIONS 


Production  of  ISAR  images  using  time-frequencies  methods  instead  of  ad¬ 
vanced  algorithms  for  motion  compensation  has  been  analyzed  in  this  thesis  because 
of  their  potential  for  fast  processing  speeds.  In  military  applications,  time  is  always 
an  important  variable.  There  are  some  situations  in  which  seconds  could  mean  lives. 
The  use  of  Fourier  and  wavelet  transforms  applied  to  ISAR  imaging  can  be  very  im¬ 
portant  for  identifying  flying  targets  in  a  very  fast  way. 

We  compared  the  ISAR  images  produced  using  the  STFT  and  the  CWT.  To 
accomplish  this  objective,  complex  backscatter  data  from  a  backhoe  loader  was  used 
to  exercise  the  algorithms. 

The  big  problem  with  the  Fourier  approach  is  its  limitation  concerning  the 
uncertainty  principle:  it’s  not  possible  to  resolve  both  time  and  frequency  with  arbi¬ 
trary  accuracy.  On  the  other  hand,  the  window  width  can  be  adjusted  for  different 
frequencies  in  the  CWT,  so  a  good  improvement  in  resolution  can  be  obtained.  With 
the  ISAR  imaging  model,  if  the  target  drastically  changes  it’s  flight  path  the  image 
will  always  display  some  artifacts  due  to  this  change.  In  this  situation,  we  observed 
that  we  can  have  some  improvement  using  the  CWT  instead  of  the  STFT.  Because 
of  its  higher  resolution  for  higher  frequencies  some  of  the  artifacts  are  eliminated, 
a  property  that  may  allow  us  to  increase  the  probability  of  identification  (observe 
Figure  40  and  Figure  37). 

If  a  target  presents  a  complex  path,  with  fast  maneuvers  the  images  will  always  be 
displayed  as  blurred  and  with  some  extra  scattering  points  due  to  aspect  change  of 
the  target.  But,  using  the  CWT  we  can  obtain  better  results  and  increase  the  quality 
of  the  image.  In  all  the  images  we  could  observed  that  the  CWT  had  superior  reso¬ 
lution  in  the  test  cases.  Since  we  can  adjust  the  scale  in  the  wavelet  transform  ,  we 
can  always  tune  the  scale  to  obtain  better  results  in  different  flight  paths,  different 
speeds  and  different  distances  from  the  radar.  We  have  more  freedom  in  this  method 
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than  in  the  Fourier  Transform  approach.  Thinking  about  the  future,  using  some  post 
processing  technics  we  can  improve  the  image.  We  can  use  noise  filters,  erosion  tech¬ 
nics  and  others.  We  should  also  exercise  the  algorithms  with  real  data.  Following 
that,  improving  the  algorithms  to  be  used  receiving  real  time  data  from  the  radar. 
This  will  allow  us  to  obtain  different  images  each  2  to  3  sec  and  with  this  data  we 
can  better  analyze  and  identify  the  target.  This  also  allow  us  to  reconstruct  a  movie 
from  the  data  to  see  each  aspect  of  the  target  during  the  time  we  were  collecting  the 
data. 
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APPENDIX  A.  MATLAB  CODE  FOR  TEST 

CASE  1 


clear  all; 

A95=getbackhoe (95 , 42) ; 
A96=getbackhoe (96 , 42) ; 
A97=getbackhoe (97 , 42) ; 
A98=getbackhoe (98 , 42) ; 
A99=getbackhoe (99 , 42) ; 
A100=getbackhoe(100,42) ; 
A101=getbackhoe (101 , 42) ; 
A102=getbackhoe(102,42) ; 

°/„elv=42.5;  94.5<Az>95.5 
B95=A95.hh(: ,1, :) ; 

Xk=  []  ; 

MN95=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B95(f ,1,1) ; 

MN95= [MN95,Xk] ; 

end 

°/„elv=42.5;  95.5<Az>96.5 
B96=A96.hh(: ,1, :) ; 

Xk=  []  ; 

MN96=  []  ; 
for  1=1 : 14 ; 
f =170 : 214; 

Xk=B96(f ,1,1) ; 

MN96= [MN96,Xk]  ; 

end 

°/0elv=42 . 5 ;  96.5<Az>97.5 
B97=A97.hh(: ,1, :) ; 

Xk=  []  ; 

MN97=  []  ; 
for  1=1 : 14 ; 
f=170 : 214; 

Xk=B97(f ,1,1) ; 

MN97= [MN97,Xk]  ; 


45 


end 


y„elv=42.5;  97.5<Az>98.5 
B98=A98.hh(: ,1, :) ; 

Xk=  []  ; 

MN98=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B98 (f , 1 , i) ; 

MN98= [MN98,Xk] ; 

end 

°/„elv=42.5;  98.5<Az>99.5 
B99=A99 . hh( : , 1 , :  )  ; 

Xk=  []  ; 

MN99=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B99 (f , 1 , i) ; 

MN99=  [MN99,Xk] ; 

end 

°/„elv=42.5;  99 . 5<Az>100 . 5 
MN100=[]  ; 

Xk=  []  ; 

B100=A100.hh(: ,1, :) ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B100 (f , 1 , i) ; 

MN100=  [MN100 , Xk] ; 

end 

°/„elv=42.5;  100 . 5<Az>101 . 5 
B101=A101.hh(: ,1, :) ; 

Xk=  []  ; 

MN101=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B101 (f , 1 , i) ; 

MN101= [MN101 , Xk] ; 

end 

°/0elv=42 . 5 ;  101 . 5<Az>102 . 5 
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B102=A102.hh(: ,1, :) ; 

MN102=[]  ; 

Xk=  []  ; 
for  i= 1 : 14 ; 
f =170 : 214; 

Xk=B102 (f , 1 , i) ; 

MN102= [MN102 , Xk] ; 

end 

°/0  Data  MAtrix  formation 

MN2  = [MN95 ,  MN96,  MN97,  MN98,  MN99,  MN100,  MN101,  MN102] ; 

7„  Display  data  matrix 
7.  p2=abs  (MN2)  ; 

%  figure,  imshow(p2 ,  [] ) ; 

7„  xlabel  ( ’theta  -  95  to  102,  phi=42deg’); 

%  ylabel (’ frequency  (amplitude)’); 

7„  Range  profiles  calculation 
DR=  []  ; 

for  i=l : 112; 

XR=f f tshif t (f f t (MN2 ( : , i) ) ) ; 

DR= [DR , XR] ; 

end 

°/0Display  range  profiles 
7,  DRABS=abs  (DR)  ; 

7o  figure,  imshow (DRABS ,[])  ; 

7o  Apply  Short-time  Fourier  transform  to  generate  ISAR  Image 
Sf=[]  ; 

L=32 ; 

noverlap=16 ; 
for  i= 1 : 45 ; 

S=spectrogram(DR(i , : ) ,hamming(80) ,40) ; 

Sf (: , : , i)=fftshift (S) ; 
end 

SFABS=abs (Sf ) ; 

70Display  ISAR  Image 

Image  =  [SFABS (:,:,!)  SFABS(:,:,2)  SFABS(:,:,3)  SFABS(:,:,4)  SFABS(:,:,5) 
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SFABS ( : , : 

,39) 

SFABS ( : , : 

,40) 

SFABS ( : , 

,41) 

SFABS ( 

3 

,42) 

SFABS ( 

, : ,43) 
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Imagew=  Image/ (max (max (Image) ) ) ; 
figure  ,  imshow(Imagew,  [0  1]); 
°/0colormap(  1-gray)  ; 
ylabel ( ’ Cross-range ’ ) ; 
xlabel ( ’ Down-Range ’ ) ; 


/(Choosing  frame  to  display 
Image2=  []  ; 
for  i=l:45; 

Imf=SFABS(: ,l,i) ; 

Image2= [Image2 , Imf ] ; 
end 


0/0Di splay  ISAR  IMAGE 
figure  ,  imshow(Image2 ,  [] ) ; 
ylabel ( ’ Cross-range ’ ) ; 
xlabel ( ’ Down-range 1 ’ ) ; 


y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

clear  all; 

A95=getbackhoe (95 , 42) 
A96=getbackhoe (96 , 42) 
A97=getbackhoe (97 , 42) 
A98=getbackhoe (98 , 42) 
A99=getbackhoe (99 , 42) 
A100=getbackhoe(100,42) ; 
A101=getbackhoe (101 , 42) ; 
A102=getbackhoe(102,42) ; 


Algorithm  2 


y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y  y 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /  0  /o  /o  /o  /o  /o  /o  /o 


7„elv=42.5;  94.5<Az>95.5 
B95=A95.hh(: ,1, :) ; 

Xk=  []  ; 

MN95=  []  ; 
for  i= 1 : 14 ; 
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f =170 : 214; 

Xk=B95 (f , 1 , i) ; 

MN95= [MN95,Xk] ; 

end 

7„elv=42.5;  95.5<Az>96.5 
B96=A96.hh(: ,1, :) ; 

Xk=  []  ; 

MN96=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B96 (f , 1 , i) ; 

MN96= [MN96,Xk] ; 

end 

7„elv=42.5;  96.5<Az>97.5 
B97=A97 . hh( : , 1 , :  )  ; 

Xk=  []  ; 

MN97=  []  ; 
for  i= 1 : 14 ; 
f =170 : 214; 

Xk=B97 (f , 1 , i) ; 

MN97= [MN97,Xk] ; 

end 

%elv=42.5;  97.5<Az>98.5 
B98=A98.hh(: ,1, :) ; 

Xk=  []  ; 

MN98=  []  ; 
for  i= 1 : 14 ; 
f =170 : 214; 

Xk=B98 (f , 1 , i) ; 

MN98= [MN98,Xk] ; 

end 

7»elv=42 . 5 ;  98.5<Az>99.5 
B99=A99 . hh( : , 1 , :  )  ; 

Xk=  []  ; 

MN99=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B99 (f , 1 , i) ; 
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MN99=  [MN99,Xk] ; 


end 

7„elv=42.5;  99 . 5<Az>100 . 5 
MN100=[]  ; 

Xk=  []  ; 

B100=A100.hh(: ,1, :) ; 
for  i=l:14; 
f =170 : 214; 

Xk=B100 (f , 1 , i) ; 

MN100=  [MN100 , Xk] ; 

end 

°/0elv=42 . 5 ;  100 . 5<Az>101 . 5 
B101=A101.hh(: ,1, :) ; 

Xk=  []  ; 

MN101=  []  ; 
for  i=l:14; 
f =170 : 214; 

Xk=B101 (f , 1 , i) ; 

MN101= [MN101 , Xk] ; 

end 

7„elv=42.5;  101 . 5<Az>102 . 5 
B102=A102.hh(: ,1, :) ; 

MN102=[]  ; 

Xk=  []  ; 
for  i=l:14; 
f =170 : 214; 

Xk=B102 (f , 1 , i)  ; 

MN102=  [MN102 , Xk] ; 

end 

%  Data  MAtrix  formation 

MN2  =  [MN95 ,  MN96,  MN97,  MN98,  MN99,  MN100,  MN101,  MN102] ; 
°/0  Display  data  matrix 
7«p2=abs  (MN2)  ; 

7.  figure,  imshow(p2 ,  []  )  ; 

7o  xlabel( ’theta  -  95  to  102,  phi=42deg’); 

7.  ylabel (’ frequency  (amplitude)’); 

7o  Range  profiles  calculation 
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RP=  []  ; 

for  i=l : 112 ; 

XR=f ft shift (f f t2(MN2 ( : , i) ) ) ; 

RP= [RP , XR] ; 

end 

°/0Pading  array 

RFP  =  padarray(RP,  [8  8]); 

Apply  Short-time  Wavelet  transform  to  generate  ISAR  Image 
Sf=[]  ; 

F=[] ; 

L=32 ; 

noverlap=16 ; 
cd=  []  ; 

for  i  =  1:61; 

°/0c  =  cwt  (RFP(i ,  :  )  ,  1 :25,  ’  sym4’ )  ; 

%  c  =  cwt  (RFP  (i, :) ,  [3  18  12.9  7  1.5],,db2’); 

%  c  =  cut  (RFP  (i,  :)  ,  [16  8  6  : -1 : 1]  ,  'morl  ’  )  ;  °/„16  8  6 
c  =  cwt  (RFP  (i,  :)  ,1:5,  ’mexh’)  ;  °/„2:5 

Sf(:,:,i)  =if  f  tshift  (if  f  t2  (c) )  ;  °/0coeff  icients 

end 

SFABS=abs (Sf ) ; 

°/oDisplay  ISAR  Image 

Image  =  [SFABS ( : , : , 1) ;  SFABS ( : , : , 2) ;  SFABS ( : , : , 3) ;  SFABS( : , : ,4) ;  SFABS(:,:,5) 
SFABS ( : , : , 6) ;  SFABS ( : , : , 7) ;  SFABS ( : , : , 8) ;  SFABS ( : , : , 9) ;  SFABS (:,:, 10) ; 

SFABS ( : , : ,15) ; 
SFABS ( : , : ,20) ; 
SFABS ( : , : ,25) ; 
SFABS ( : , : ,30) ; 
SFABS ( : , : ,35) ; 
SFABS ( : , : ,40) ; 
SFABS ( : , : ,45)] ; 

Imagew=  Image/ (max (max (Image) )) ; 

°/0Imag  =  imadjust (Imagew,  [.07;  .5]  ,  [0;  1] )  ; 

°/0Imag  =  histeq(Imagew, 64)  ; 

°/0Imag  =  adapthisteq(Imageu)  ; 
figure,  imshow( Imagew,  [] ) ; 

%  colormap ( 1-gray) ; 
ylabel ( ’ Range ’ ) ; 
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xlabel ( ’ Down-range ’ ) ; 

°/0Choosing  frame  to  display 
Image2=  []  ; 
for  i=l:45; 

Imf 2=SFABS (2 , : ,i) ; 

Image2= [Image2 ; Imf 2] ; 
end 

“/.Display  ISAR  IMAGE 
figure  ,  imshow(Image2 , [] ) 
ylabel ( ’ Cross-range ’ ) ; 
xlabel ( ’ Down-range ’ ) ; 


APPENDIX  B.  MATLAB  CODE  FOR  TEST 

CASE  2 


clear  all; 

7„  Extraction  data  from  database 
A95=getbackhoe (95 , 38) ; 
A96=getbackhoe (96 , 38) ; 
A97=getbackhoe (97 , 38) ; 
A98=getbackhoe (98 , 38) ; 
A99=getbackhoe (99 , 39) ; 
A100=getbackhoe (99 , 40) ; 
A101=getbackhoe (99 , 41) ; 
A102=getbackhoe (99 , 42) ; 

°/„elv=38.5;  94.5<Az>95.5 
B95=A95.hh(: ,1, :) ; 

Xk=  []  ; 

MN95=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B95 (f , 1 , i) ; 

MN95= [MN95,Xk] ; 

end 

°/„elv=38.5;  95.5<Az>96.5 
B96=A96.hh(: ,1, :) ; 

Xk=  []  ; 

MN96=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B96 (f , 1 , i) ; 

MN96= [MN96,Xk] ; 

end 

°/„elv=38 . 5 ;  96 . 5<Az>97 . 5 
B97=A97.hh(: ,1, :) ; 

Xk=  []  ; 

MN97=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B97 (f , 1 , i) ; 
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MN97=  [MN97,Xk] ; 


end 

°/„elv=38 . 5 ;  97 . 5<Az>98 . 5 
B98=A98 . hh( : , 1 , :  )  ; 

Xk=  []  ; 

MN98=  []  ; 
for  i=l:14; 
f =170 : 214; 

Xk=B98 (f , 1 , i) ; 

MN98= [MN98,Xk] ; 

end 

°/038 . 5<elv>39 . 5 ;  Az=99.5 
B99=A99.hh(: , : ,7) ; 

Xk=  []  ; 

MN99=  []  ; 
for  i=l:14; 
f =170 : 214; 

Xk=B99(f ,i,l) ; 

MN99= [MN99,Xk] ; 

end 

°/039 . 5<elv>40 . 5 ;  Az=99.5 

MN100=[]  ; 

Xk=  []  ; 

B100=A100.hh(: , : ,7) ; 
for  i=l:14; 
f =170 : 214; 

Xk=B100 (f , i ,  1)  ; 
MN100=  [MN100 , Xk] ; 

end 

%40 . 5<elv>41 . 5 ;  Az=99.5 
B101=A101.hh(: , : ,7) ; 

Xk=  []  ; 

MN101=  []  ; 
for  i= 1 : 14 ; 
f =170 : 214; 

Xk=B101 (f , i , 1) ; 
MN101= [MN101 , Xk] ; 
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end 


°/„41 . 5<elv>42 . 5 ;  Az=99.5 

B102=A102 ,hh( : , : ,7) ; 

MN102=[]  ; 

Xk=  []  ; 
for  i= 1 : 14 ; 
f =170 : 214; 

Xk=B102 (f , i , 1) ; 

MN102= [MN102 , Xk] ; 

end 

°/0  Data  Matrix  formation 

MN2  = [MN95 ,  MN96,  MN97,  MN98,  MN99,  MN100,  MN101,  MN102] ; 

°/0Plot  data  Matrix 

p2=abs (MN2) ; 

figure,  imshow(p2,  [] ) ; 

xlabel( ’theta  -  95  to  102,  phi=42deg’); 

ylabel (’ frequency  (amplitude)’); 

°/0Range  profiles  calculation 
DR=  []  ; 

for  i=l : 112; 

XR=f f tshif t (f f t (MN2 ( : , i) ) ) ; 

DR= [DR , XR] ; 

end 

7„  Display  Range-profiles 
°/„  DRABS=abs  (DR)  ; 

%  figure,  imshow  (DRABS ,[])  ; 

%  Apply  Fourrier  tranform  to  all  lines  to  form  the  ISAR  image 
Sf=[]  ; 

L=32 ; 

noverlap=16 ; 
for  i= 1 : 45 ; 

S=spectrogram(DR(i , : ) , hamming (100) , 50) ; 

Sf (: , : , i)=fftshift (S) ; 
end 

SFABS=abs (Sf ) ; 
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Image  =  [SFABS ( : , : , 1)  SFABS(:,:,2)  SFABS(:,:,3)  SFABS(:,:,4)  SFABS(:,:,5) 
SFABS ( : , : , 6)  SFABS(:,:,7)  SFABS(:,:,8)  SFABS(:,:,9)  SFABS 10) 
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Imagew=  Image/ (max (max (Image) )) ; 
figure,  imshow(Imagew, [0  1]); 
ylabel (’ Cross  Range’); 
xlabel ( ’ Range ’ )  ; 

7„  Choose  one  frame  to  display 
Image2=  []  ; 
for  i= 1 : 45 ; 

Imf =SFABS ( : , 1 , i) ; 

Image2= [Image2 , Imf ] ; 
end 


7„  Display  frame 
figure,  imshow(Image2 ,  [] ) ; 
ylabel ( ’ Range ’ ) ; 
xlabel (’ Cross  Range’); 


0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0/  0. 


Algorithm  2  nnunnununiinun 


clear  all; 

A95=getbackhoe (95 , 38) ; 
A96=getbackhoe (96 , 38) ; 
A97=getbackhoe (97 , 38) ; 
A98=getbackhoe (98 , 38) ; 
A99=getbackhoe (99 , 39) ; 
A100=getbackhoe (99 , 40) ; 
A101=getbackhoe (99 , 41) ; 
A102=getbackhoe (99 , 42) ; 

°/„elv=38.5;  94.5<Az>95.5 
B95=A95.hh(: ,1, :) ; 
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Xk=  []  ; 

MN95=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B95 (f , 1 , i) ; 

MN95= [MN95,Xk] ; 

end 

7„elv=38.5;  95.5<Az>96.5 
B96=A96.hh(: ,1, :) ; 

Xk=  []  ; 

MN96=  []  ; 
for  i=l : 14; 
f=170 : 214; 

Xk=B96 (f , 1 , i) ; 

MN96= [MN96,Xk] ; 

end 

7„elv=38 . 5 ;  96 . 5<Az>97 . 5 
B97=A97 . hh( : , 1 , :  )  ; 

Xk=  []  ; 

MN97=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B97 (f , 1 , i) ; 

MN97= [MN97,Xk] ; 

end 

7.elv=38 . 5 ;  97 . 5<Az>98 . 5 
B98=A98.hh(: ,1, :) ; 

Xk=  []  ; 

MN98=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B98 (f , 1 , i) ; 

MN98= [MN98,Xk] ; 

end 

7«38 . 5<elv>39 . 5 ;  Az=99.5 
B99=A99.hh(: , : ,7) ; 

Xk=  []  ; 

MN99=  []  ; 
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for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B99 (f  ,  i , 1) ; 

MN99=  [MN99,Xk] ; 

end 

°/039 . 5<elv>40 . 5 ;  Az=99.5 
MN100=[]  ; 

Xk=  []  ; 

B100=A100.hh(: , : ,7) ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B100 (f , i , 1) ; 

MN100=  [MN100 , Xk] ; 

end 

7040 . 5<elv>41 . 5 ;  Az=99.5 
B101=A101.hh(: , : ,7) ; 

Xk=  []  ; 

MN101=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B101 (f , i , 1) ; 

MN101=  [MN101 , Xk] ; 

end 

°/„41 . 5<elv>42 . 5 ;  Az=99.5 
B102=A102 . hh( : , : ,7) ; 

MN102=  []  ; 

Xk=  []  ; 
for  i= 1 : 14 ; 
f=170 : 214; 

Xk=B102 (f , i , 1) ; 

MN102= [MN102 , Xk] ; 

end 

%  Data  Matrix  formation 

MN2  =  [MN95 ,  MN96,  MN97,  MN98,  MN99,  MN100,  MN101,  MN102] 

7«Di splay  data  matrix 

p2=abs (MN2) ; 

figure,  imshow(p2,  [] ) ; 
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xlabel( ’theta  -  95  to  102,  phi=42deg’); 
ylabel (’ frequency  (amplitude)’); 

°/0Range  profiles  calculation 
RP=  []  ; 

for  i=l : 112 ; 

XR=f ft shift (f f t2(MN2 ( : , i) ) ) ; 

RP= [RP , XR] ; 

end 

°/0  Array  Padding 

RFP  =  padarray(RP,  [8  8]); 

°/o7o  Apply  Wavelet  tranform  to  all  lines  to  form  the  ISAR  image 
Sf=[]  ; 

F=[]  ; 

L=32 ; 

noverlap=16 ; 
cd=  []  ; 

for  i  =  1:61; 

°/0c  =  cwt  (RFP(i ,  : )  ,  1 :  25 ,  ’  sym4’ )  ; 
c  =  cwt (RFP (i , : ) ,  [3  18  12.9  7  1.5],’db2’); 

7,  c  =  cwt  (RFP  (i,  :)  ,  [16  8  6  : -1 : 1]  , ’morl  ’  )  ; 

7.  c  =  cwt  (RFP  (i ,  : )  ,  10  : 13 , ’mexh’ )  ; 

Sf(:,:,i)  =if  f  tshif  t  (if  ft2  (c) )  ;  7«coef  f  icients 
end 

SFABS=abs (Sf ) ; 

7oDisplay  ISAR  image 

Image  =  [SFABS ( : , : , 1) ;  SFABS ( : , : , 2) ;  SFABS ( : , : , 3) ;  SFABS( : , : ,4) ;  SFABS ( : , : , 5) ; 
SFABS ( : , : , 6) ;  SFABS( : , : ,7) ;  SFABS ( : , : , 8) ;  SFABS ( : , : , 9) ;  SFABS (:,:, 10) ; 
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Imagew=  Image/ (max (max (Image) )) ; 
figure,  imshow(Imagew, [] ) ; 
ylabel ( ’ Range 1 ’ )  ; 
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xlabel(’ Doppler  frequency’); 


“/Display  ISAR  Image 
Image2=  []  ; 
for  i=l : 45 ; 
Imf2=SFABS(l, : ,i) ; 
Image2= [Image2 ; Imf 2] ; 
end 
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